
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
#Map of Europe
europe_map =plt.imread('map.png')

#Lists of city coordinates
city_coords={"Barcelona":[2.154007, 41.390205], "Belgrade": [20.46,44.79], "Berlin": [13.40,52.52], "Brussels":[4.35,50.85],"Bucharest":[26.10,44.44], "Budapest": [19.04,47.50], "Copenhagen":[12.57,55.68], "Dublin":[-6.27,53.35], "Hamburg": [9.99, 53.55], "Istanbul": [28.98, 41.02], "Kiev": [30.52,50.45], "London": [-0.12,51.51], "Madrid": [-3.70,40.42], "Milan":[9.19,45.46], "Moscow": [37.62,55.75], "Munich": [11.58,48.14], "Paris":[2.35,48.86], "Prague":[14.42,50.07], "Rome": [12.50,41.90], "Saint Petersburg": [30.31,59.94], "Sofia":[23.32,42.70], "Stockholm": [18.06,60.33],"Vienna":[16.36,48.21],"Warsaw":[21.02,52.24]}


#Helper code for plotting plans
#First, visualizing the cities.
import csv
with open("european_cities.csv", "r") as f:
    data = list(csv.reader(f, delimiter=';'))
    cities = data[0]
    
fig, ax = plt.subplots(figsize=(10,10))

ax.imshow(europe_map, extent=[-14.56,38.43, 37.697 +0.3 , 64.344 +2.0], aspect = "auto")



# Map (long, lat) to (x, y) for plotting
for city,location in city_coords.items():
    x, y = (location[0], location[1])
    plt.plot(x, y, 'ok', markersize=5)
    plt.text(x, y, city, fontsize=12);


#A method you can use to plot your plan on the map.
def plot_plan(city_order):
    fig, ax = plt.subplots(figsize=(10,10))
    ax.imshow(europe_map, extent=[-14.56,38.43, 37.697 +0.3 , 64.344 +2.0], aspect = "auto")

    # Map (long, lat) to (x, y) for plotting
    for index in range(len(city_order) -1):
        current_city_coords = city_coords[city_order[index]]
        next_city_coords = city_coords[city_order[index+1]]
        x, y = current_city_coords[0], current_city_coords[1]
        #Plotting a line to the next city
        next_x, next_y = next_city_coords[0], next_city_coords[1]
        plt.plot([x,next_x], [y,next_y])
        
        plt.plot(x, y, 'ok', markersize=5)
        plt.text(x, y, index, fontsize=12);
    #Finally, plotting from last to first city
    first_city_coords = city_coords[city_order[0]]
    first_x, first_y = first_city_coords[0], first_city_coords[1]
    plt.plot([next_x,first_x],[next_y,first_y])
    #Plotting a marker and index for the final city
    plt.plot(next_x, next_y, 'ok', markersize=5)
    plt.text(next_x, next_y, index+1, fontsize=12);
    plt.show();


#Example usage of the plotting-method.
plan = list(city_coords.keys()) # Gives us the cities in alphabetic order
print(plan)
plot_plan(plan)

['Barcelona', 'Belgrade', 'Berlin', 'Brussels', 'Bucharest', 'Budapest', 'Copenhagen', 'Dublin', 'Hamburg', 'Istanbul', 'Kiev', 'London', 'Madrid', 'Milan', 'Moscow', 'Munich', 'Paris', 'Prague', 'Rome', 'Saint Petersburg', 'Sofia', 'Stockholm', 'Vienna', 'Warsaw']


# Implement the algorithm here
import itertools

def exhaustiveSearch(n):
    # Create a dictionary of distances between cities
    distances = {}
    for i in range(n):
        for j in range (n):
            # Only calculate distance between different cities
            if cities[i] != cities[j]:
                distances[(cities[i], cities[j])] = data[i+1][j]

    # Get a list of all permutations of n cities
    nCities = cities[:n]
    combs = list(itertools.permutations(nCities))

    # Set initial values for shortest distance and shortest path
    shortest = float("inf")
    shortestPath = None

    # Loop through each possible path and calculate the total distance    shortest = float("inf")
    shortestPath = None
    for path in combs:
        dist = 0
        for i in range(n-1):
            city = path[i]
            city2 = path[i+1]
            dist += float(distances[(city, city2)])
            if dist > shortest:
                # Stop calculating the distance if it's already longer than the current shortest
                break
        # Add the distance from the last city to the first city to complete the path
        dist += float(distances[(city2, path[0])])
        # Update the shortest path and distance if this is shorter
        if dist < shortest:
            shortest = dist
            shortestPath = path
    return shortestPath, shortest


# Answer
import time

start = time.time()
path, distance = exhaustiveSearch(10)
end = time.time()
ES_time = end - start
print("Executed time for Exhaustive search: ",  "{:.4f}".format(ES_time))

print("The shortest tour found : ", path, "with distance ", "{:.4f}".format(distance))
plot_plan(path)

Executed time for Exhaustive search:  16.9549
The shortest tour found :  ('Copenhagen', 'Hamburg', 'Brussels', 'Dublin', 'Barcelona', 'Belgrade', 'Istanbul', 'Bucharest', 'Budapest', 'Berlin') with distance  7486.3100


# Implement the algorithm here
import random
# Collecting distances in dict
distances = {}
for i in range(24):
    for j in range (24):
        if cities[i] != cities[j]:
            distances[(cities[i], cities[j])] = data[i+1][j]

# Method for calculating dist (also used in part3)
def get_distance(path, n=24, best = float("inf")):
        dist = 0
        for i in range(n-1):
            city = path[i]
            city2 = path[i+1]
            dist += float(distances[(city, city2)])
            if dist > best:
                return              # quits when path is worse
        dist += float(distances[(city2, path[0])])
        return dist

# This function returns the best tour and its distance 
# for a given number of cities n
def hill_climbing(n):
    nCities = cities[:n]        # global list
    
    # Swap 2 neighboring cities in given neighbor path
    def find_neighbor(nb):
        i = random.randint(0,n-2)
        j = i + 1
        nb[i], nb[j] = nb[j], nb[i]
        return nb

    best_tour = nCities
    random.shuffle(best_tour)
    best = get_distance(best_tour, n)
    nb_path = best_tour  
    visited_neighbors = []

    i = 0
    while(i < 1000000):  #avoid infinite loops
        nb = find_neighbor(nb_path)
        # Avoid the visited neighbors
        if not nb_path in visited_neighbors:
            visited_neighbors.append(nb_path[:])    #adding copy
            nb_dist = get_distance(nb_path, n, best)

            #  Update best path if better
            if nb_dist != None:     
                best = nb_dist
                best_tour = nb_path
            else:
                break           #if new path is not better
        i += 1
    
    return best_tour, best

print("HILL-CLIMBING FOR 10 CITIES")
start = time.time()
all_tours = {}
for i in range(20):
    tour, tour_d = hill_climbing(10)
    all_tours[tour_d] = tour
end = time.time()
HC_time = end - start

print("Executed time for 20 runs of the hill climbing method: ", "{:.10f}".format(HC_time))
print("Compared to Exhaustive Search, this method performed" , "\033[1m{:.4f}\033[0m".format(ES_time - HC_time), "s better!")


keys = all_tours.keys()
best, worst, avg, std = min(keys), max(keys), sum(keys)/len(keys), np.std(list(keys))
best_tour = all_tours.get(best)
print("\nThe plot shows the best tour of the 20 runs: ", "{:.4f}".format(best), "\nWith path : ", best_tour)
plot_plan(best_tour)

print("The worst path of the 20 runs had a distance of : ", "{:.4f}".format(worst))
print("The mean: ", "{:.4f}".format(avg))
print("The standard deviation: ", "{:.4f}".format(std))

print("\nHILL-CLIMBING FOR 24 CITIES")
start = time.time()
all_tours = {}
for i in range(20):
    tour, tour_d = hill_climbing(24)
    all_tours[tour_d] = tour
end = time.time()
HC_time = end - start

print("Executed time for 20 runs of the hill climbing method: ", "{:.10f}".format(HC_time))
print("Compared to Exhaustive Search, this method performed" , "\033[1m{:.4f}\033[0m".format(ES_time - HC_time), "s better!")


keys = all_tours.keys()
best, worst, avg, std = min(keys), max(keys), sum(keys)/len(keys), np.std(list(keys))
best_tour = all_tours.get(best)
print("\nThe plot shows the best tour of the 20 runs: ", "{:.4f}".format(best), "\nWith path : ", best_tour)
plot_plan(best_tour)

print("The worst path of the 20 runs had a distance of : ", "{:.4f}".format(worst))
print("The mean: ", "{:.4f}".format(avg))
print("The standard deviation: ", "{:.4f}".format(std))

HILL-CLIMBING FOR 10 CITIES
Executed time for 20 runs of the hill climbing method:  0.0009448528
Compared to Exhaustive Search, this method performed 16.9540 s better!

The plot shows the best tour of the 20 runs:  9045.7700 
With path :  ['Belgrade', 'Budapest', 'Barcelona', 'Hamburg', 'Brussels', 'Dublin', 'Bucharest', 'Berlin', 'Copenhagen', 'Istanbul']

The worst path of the 20 runs had a distance of :  13940.7300
The mean:  11374.1070
The standard deviation:  1272.9970

HILL-CLIMBING FOR 24 CITIES
Executed time for 20 runs of the hill climbing method:  0.0000000000
Compared to Exhaustive Search, this method performed 16.9549 s better!

The plot shows the best tour of the 20 runs:  27060.7500 
With path :  ['Belgrade', 'London', 'Milan', 'Stockholm', 'Paris', 'Vienna', 'Saint Petersburg', 'Budapest', 'Istanbul', 'Barcelona', 'Bucharest', 'Madrid', 'Dublin', 'Copenhagen', 'Munich', 'Brussels', 'Prague', 'Berlin', 'Hamburg', 'Kiev', 'Rome', 'Warsaw', 'Moscow', 'Sofia']

The worst path of the 20 runs had a distance of :  33851.6500
The mean:  30378.8270
The standard deviation:  2035.5661


# Implement the algorithm here

from collections import OrderedDict

def pmx_pair(p1, p2):
    l = len(p1)
    # choose segment
    left = random.randint(0,len(p1))
    right = random.randint(0,len(p1))
    if left > right:
        left, right = right, left

    # Initialize children with only the segments
    c1 = [None] * l
    c1[left:right] = p2[left:right]
    c2 = [None] * l
    c2[left:right] = p1[left:right]

    # Loop through parent 1
    for i in range(len(p1)):
        if not c1[i]:      # if the value is not yet assigned to the child 1
            j = p1[i]       
            idx = None
            while j in c1:         # if the value is already in the child 1
                idx = c1.index(j)
                j = c2[idx]         # find the corresponding value from child 2
            c1[i] = j

    # Doing the same for parent 2 with child 2
    for i in range(len(p2)):    
        if not c2[i]:      
            j = p2[i] 
            idx = None
            while j in c2:         
                idx = c2.index(j)
                j = c1[idx]
            c2[i] = j
            
    return c1, c2

# Method for finding children with pmx
# returns a list of parents and children
def find_offsprings(sols):
    for i in range(0, len(sols), 2):
        child1, child2 = pmx_pair(sols[i], sols[i+1])
        if child1 not in sols:
            sols.append(child1)
        if child2 not in sols:
            sols.append(child2)
    return sols


def genetic(n, p, MAX, plotY):
    solutions = []
    lst = cities[:n].copy()
    # Initialization of n parents by random shuffling
    for i in range(p):         
        random.shuffle(lst)
        solutions.append(lst[:])

    # Starting recombination for given generations
    gen = 0
    while (gen < MAX):
        solutions = find_offsprings(solutions)

        shortest_paths = sorted(solutions, key=lambda i: get_distance(i, n))   # sorts by distance
        solutions = shortest_paths[:p]      # keep n best from this gen
        
        best_path = solutions[0]
        best = get_distance(best_path, n)

        plotY[gen] += best         # adds the best distance to the plot data
        gen += 1

    worst = get_distance(solutions[p-1], n)
    return plotY, (best, best_path), worst      # saves the best solution in tuple

# Running EA on n parents for certain generations
def runEA(n, p, gens):
    start = time.time()
    plotY = [0] * gens
    best_sols = []
    best, worst = (float("inf"), None), 0

    # Run 20 times
    for i in range(20):
        plotY, b, w = genetic(n, p, gens, plotY)
        # Find the best & worst
        if min(best[0], b[0]) == b[0]:
            best = b
        best_sols.append(best[0])
        worst = max(worst, w)

    # Finds the average fitness over 20 runs
    plotY = np.array(plotY)/20
    end = time.time()

    print("\nThe evolutionary algorithm for population size", p, "was performed in ", round(end - start, 3), "s.")
    print("The best solution gives: ", round((best[0]), 3))
    print("The worst solution gives: ", round((worst), 3))    
    print("The mean of the best solutions gives: ", round((np.mean(best_sols)), 3))
    print("The standard deviation of the best solutions gives: ", round((np.std(best_sols)), 3))
    return plotY, best[1]      #returns the datapoints and the best path

gens = 100
pops = [100, 300, 500]
plt.title("Average fitness over runs for every gen")
plotX = np.linspace(0, gens, gens)      # x-axis

# Running EA for 24 cities on every population size
for i in range(len(pops)):
    plotY, best = runEA(24, pops[i], gens)
    pops[i] = best
    plt.plot(plotX, plotY)

# Showing plots of best solution of every population
for i in range(len(pops)):
    plot_plan(pops[i])
    print("This is the best solution for ", i+1, ". test with path: ", pops[i])

print("Conclusion: The best of them all is the test with size 300, because the best solution isn't much different from size 500.")
print("Even though it finds a slightly better solution, it is at the expense of performing time.")
print("At an attempt on size 600 and additional 10s, it also didn't yield a significant improvement")

The evolutionary algorithm for population size 100 was performed in  5.459 s.
The best solution gives:  12639.41
The worst solution gives:  16630.79
The mean of the best solutions gives:  13138.05
The standard deviation of the best solutions gives:  521.372

The evolutionary algorithm for population size 300 was performed in  18.761 s.
The best solution gives:  12504.65
The worst solution gives:  15140.46
The mean of the best solutions gives:  12575.126
The standard deviation of the best solutions gives:  56.328

The evolutionary algorithm for population size 500 was performed in  38.445 s.
The best solution gives:  12325.93
The worst solution gives:  14929.04
The mean of the best solutions gives:  12421.297
The standard deviation of the best solutions gives:  118.215

This is the best solution for  1 . test with path:  ['Paris', 'Brussels', 'Prague', 'Berlin', 'Hamburg', 'Copenhagen', 'Stockholm', 'Saint Petersburg', 'Moscow', 'Kiev', 'Bucharest', 'Istanbul', 'Sofia', 'Belgrade', 'Budapest', 'Warsaw', 'Vienna', 'Munich', 'Milan', 'Rome', 'Barcelona', 'Madrid', 'Dublin', 'London']

This is the best solution for  2 . test with path:  ['Saint Petersburg', 'Stockholm', 'Copenhagen', 'Hamburg', 'Brussels', 'Paris', 'London', 'Dublin', 'Madrid', 'Barcelona', 'Milan', 'Rome', 'Munich', 'Prague', 'Berlin', 'Warsaw', 'Vienna', 'Budapest', 'Belgrade', 'Sofia', 'Istanbul', 'Bucharest', 'Kiev', 'Moscow']

This is the best solution for  3 . test with path:  ['Sofia', 'Belgrade', 'Budapest', 'Vienna', 'Warsaw', 'Prague', 'Munich', 'Milan', 'Rome', 'Barcelona', 'Madrid', 'Dublin', 'London', 'Paris', 'Brussels', 'Hamburg', 'Berlin', 'Copenhagen', 'Stockholm', 'Saint Petersburg', 'Moscow', 'Kiev', 'Bucharest', 'Istanbul']
Conclusion: The best of them all is the test with size 300, because the best solution isn't much different from size 500.
Even though it finds a slightly better solution, it is at the expense of performing time.
At an attempt on size 600 and additional 10s, it also didn't yield a significant improvement


# Answer
# Running EA for 10 cities on every population size
pops = [100, 300, 500]
for i in range(len(pops)):
    plotY, best = runEA(10, pops[i], gens)
    pops[i] = best
    plt.plot(plotX, plotY)

# Showing plots of best solution of every population
for i in range(len(pops)):
    plot_plan(pops[i])
    print("This is the best solution for ", i+1, ". test with path: ", pops[i])

print("My GA found the shortest tour for 10 cities with exact same distance.")
print("For 10 cities, the exhaustive search did better than mine with the population size 500,\n but size 100 was more than good enough to find the best solution.")
print("The performing time was only 3s, 13s better!")
print("For 24 cities, the result always varied for everytime I tried running it. Therefore exhaustive search would have found the best solution")
print("at first attempt, but at the expense of a lot more time. Hence GA is much more effective, but not always optimal to find the best solution.")

The evolutionary algorithm for population size 100 was performed in  4.197 s.
The best solution gives:  7486.31
The worst solution gives:  8066.12
The mean of the best solutions gives:  7486.31
The standard deviation of the best solutions gives:  0.0

The evolutionary algorithm for population size 300 was performed in  11.705 s.
The best solution gives:  7486.31
The worst solution gives:  8009.86
The mean of the best solutions gives:  7486.31
The standard deviation of the best solutions gives:  0.0

The evolutionary algorithm for population size 500 was performed in  24.331 s.
The best solution gives:  7486.31
The worst solution gives:  8150.86
The mean of the best solutions gives:  7486.31
The standard deviation of the best solutions gives:  0.0

This is the best solution for  1 . test with path:  ['Copenhagen', 'Hamburg', 'Brussels', 'Dublin', 'Barcelona', 'Belgrade', 'Istanbul', 'Bucharest', 'Budapest', 'Berlin']

This is the best solution for  2 . test with path:  ['Copenhagen', 'Hamburg', 'Brussels', 'Dublin', 'Barcelona', 'Belgrade', 'Istanbul', 'Bucharest', 'Budapest', 'Berlin']

This is the best solution for  3 . test with path:  ['Copenhagen', 'Hamburg', 'Brussels', 'Dublin', 'Barcelona', 'Belgrade', 'Istanbul', 'Bucharest', 'Budapest', 'Berlin']
My GA found the shortest tour for 10 cities with exact same distance.
For 10 cities, the exhaustive search did better than mine with the population size 500,
 but size 100 was more than good enough to find the best solution.
The performing time was only 3s, 13s better!
For 24 cities, the result always varied for everytime I tried running it. Therefore exhaustive search would have found the best solution
at first attempt, but at the expense of a lot more time. Hence GA is much more effective, but not always optimal to find the best solution.


# Implement algorithm here

	Barcelona	Belgrade	Berlin	Brussels	Bucharest	Budapest
Barcelona	0	1528.13	1497.61	1062.89	1968.42	1498.79
Belgrade	1528.13	0	999.25	1372.59	447.34	316.41
Berlin	1497.61	999.25	0	651.62	1293.40	1293.40
Brussels	1062.89	1372.59	651.62	0	1769.69	1131.52
Bucharest	1968.42	447.34	1293.40	1769.69	0	639.77
Budapest	1498.79	316.41	1293.40	1131.52	639.77	0

IN3050/IN4050 Mandatory Assignment 1 - Annica - @annicas : Traveling Salesman Problem¶

Rules¶

Delivery¶

What to deliver?¶

Introduction¶

Problem¶

Helper code for visualizing solutions¶

Exhaustive Search¶

Hill Climbing¶

Genetic Algorithm¶

Hybrid Algorithm (IN4050 only)¶

Lamarckian¶

Baldwinian¶

Task¶